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Abstract 

We investigate the phase of the quark determinant with finite chemical potential in lattice QCD 
using both analytic and numerical methods. Applying the winding number expansion and the 
hopping parameter expansion to the logarithm of the determinant, we show that the absolute value 
of the phase has an upper bound that grows with the spatial volume but decreases exponentially 
with an increase in the temporal extent of the lattice. This analytic but approximate result is 
confirmed with a numerical study in four-flavor QCD in which the phase is calculated exactly. 
Since the phase is well controlled on lattices with larger time extents, we try the phase reweighting 
method in a region beyond /u/T = 1 where the Taylor expansion method cannot be applied. 
Working in four-flavor QCD, we find a first-order like behavior on a 6^ x 4 lattice at /u/T ~ 0.8 
which was previously observed by Kentucky group with the canonical method. We also show that 
the winding number expansion has a nice convergence property beyond /i/T = 1. We expect 
that this expansion is useful to study the high density region of the QCD phase diagram at low 
temperatures. 

PACS numbers: ll.lO.Wx, 12.38.Gc 
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I. INTRODUCTION 



Establishing the QCD phase diagram spanned by the temperature T and the quark 
chemical potential in a quantitative way is an important task of lattice QCD. See Refs. [l|, 0] 
for recent progress. The Monte Carlo simulation technique, which has been successfully 
applied to the finite temperature phase transition studies in lattice QCD, cannot be directly 
applied to the finite density case due to the complexity of the quark determinant det D = 
I detD|e*^ for /x 7^ 0. A simple way out of the problem is the reweighting method which 
incorporates the absolute value of the determinant | detDl into the integral measure. This 
bipass, however, suffers from the sign problem with increasing /i caused by increasingly 
larger gauge fluctuations of the reweighting factor e*^. Understanding the property of the 
phase is crucial to control the sign problem. 

The average of the phase factor (e*^) was investigated in the framework of chiral pertur- 
bation theory in Refs. (sj-lsj and by random matrix model in Ref. jsj thus far. The phase 
itself 9 has also been investigated by Taylor expansion in lattice QCD In this article 

we report on our attempt to understand the property of the phase directly in lattice QCD. 
Near the phase transition, low energy effective theories may not be reliable. To survey such 
a region, it is imperative to work within the first principle framework of QCD itself. 

The Taylor expansion method works directly with QCD. On the T-fi plane of QCD, 
however, we expect the expansion to be reliable only in the region /i/T < 1. Therefore, 
alternative approximation methods which works for wider region in the parameter /x/T is a 
welcome step to understand the phase diagram. In this paper, we show that the winding 



number expansion whose convergence is better for larger temporal size of lattice, can 
be a possible candidate of the approximation. 

This expansion, combined with the hopping parameter expansion, allows us an analytical 
investigation of the lattice parameter dependence of the phase for the determinant of the 
Wilson operator. The result, which is obtained at the leading order of the hopping parameter 
expansion, tells us that the magnitude of the phase decreases for larger temporal size of the 
lattice with other parameters held fixed. A numerical test in A^f = 4 QCD with an exact 
calculation of the phase has confirmed that this property holds true beyond the leading 
order of the hopping parameter expansion. 

These analyses drive us to try a further numerical study with the reweighting method 
since we can avoid the sign problem by increasing the temporal size of the lattice. As a first 
testing ground, we carry out A^f = 4 QCD simulation and observe a first-order like behavior 
which was previously reported in Ref. fl2| . 

The rest of the paper is organized as follows. In Sec. [Tll we discuss the convergence 
of the winding number expansion, and the behavior of the phase as a function of lattice 
parameters. Simulation results in A^f = 4 QCD are shown in Sec. IIII CI Finally, we present 
concluding remarks. 

Throughout this paper we consider a 4-dimensional Euclidean lattice of a size A^l ^ 
A^x with the periodic (anti-periodic) boundary conditions in the time direction for boson 
(fermion) fields. A summary of the notations and definitions relevant for the winding number 
expansion and a reduction process of the quark determinant is given in appendix |Al 
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II. PHASE OF QUARK DETERMINANT 



A. Winding number expansion 



The starting point of our discussion is the winding number expansion of the logarithm 
of the quark determinant for Wilson-type fermion action [ll| . A brief review of the method 
is provided in appendix El for the case of the unimproved Wilson fermion action. Inclusion 
of the clover term is straightforward, however, and the simulation results presented later in 
this article are obtained for the clover-improved fermion action. 

For the one-flavor case, the winding number expansion takes the form. 



(1) 



where the dependence on the hopping parameter k, = l/(8 + 2amo) and the quark chemical 
potential /i is given explicitly. In the original paper the authors used T^'^^ instead of 
yi<i)_ J" = l/(^aNT) is the temperature. As explained in appendix |Al the above expression is 
derived by casting the Wilson-Dirac operator into blocks connecting time slices, and reducing 
the quark determinant in terms of these blocks. 

The first factor Aq defined in eq. flA24p is composed from mainly the block diagonal parts 
of the Wilson operator in a time blocked form. Hence it depends on k, but is independent 
of /i. This factor can be shown to be real, but is not guaranteed to be positive. On the 
other hand, V^'^\k) {q ^ 0) in the second factor, which are defined in eq. flA33p . contain 
global information in the time direction, being composed from the off-diagonal parts of the 
Wilson operator. Namely, V^'^\k) is a sum of quark loops winding around the time direction 
q times, and so appears with the g-th power of fugacity e^^'^ . It has the following properties: 



y("<?) = yiQT for 



(2) 
(3) 



The details are given in appendix \M 

Let us express the phase of the determinant in terms of the winding number expansion. 
To this end, using eqs. ([2]) and ([3]), we rewrite. 



det[D{fi)] = Aqc exp 



(2cosh(g/i/r)Re[1/('?)] +22sinh(g/i/T)Im[1/('?)]) 

9=1 



and defining the phase through 



det[D(/i)] = I det [£)(//)] I e 



we read 



det[D(/i)]| = |Ao|e-^""exp 



^2cosh(g/i/T)Re[l^('')] 



9=1 



e{^l) = arg(Ao) - ^2sinh(g/i/T)Im[\/(«)]. 

9=1 



(4) 

(5) 

(6) 
(7) 
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Since Aq is a real number, the phase contribution arg(Ao) takes the value or vr. For 
even number of degenerate flavors, it is zero because of a.Tg{Al) = 0. Even for odd number of 
flavors, we expect that arg(Ao) = holds if the corresponding quark mass is heavy enough, 
since we know that the Wilson fermion determinant with yU, = is effectively positive in 
the strange mass region. Therefore in the following we ignore the phase contribution from 
the local dynamics arg(Ao) and exclusively consider the essential part of the phase which 
depends on vanishing at /i = 0, 



e{fi) = -^2sinh(g/i/T)Im[\/('?)]. 



(8) 



q=l 



For later use, we write down the phase for the general A'^f-fiavor case. From the full 
determinant, 



l[det[D{nf,fif)] 



f 



l[det[D{Kf,fif)] 



exp 



(9) 



we read the total phase 



/ / 9=1 



(10) 



B. Convergence of the winding number expansion 

In order to identify an effective expansion parameter in the winding number expansion 
eq. dH]), we adopt two simplifications. The first factor behaves as 

2sinh(gyu/T) ~ exp(gyu/T), (11) 

for large q. For the second term, one expects 

V^"^ oc (2k)«^t (12) 

for small k and large q, since the quark loops in V^''^ have lengths of at least gA^x and each 
hop along the loop comes together with a factor of 2k,. Writing the proportionality constant 
in eq. f lT2|) as Cq, the phase is estimated as 

oo oo 

e ^ -^e('"^/^) ■c,(2«:)''^- = _^c,{e'^/^+^-^'^(2'^)}''. (13) 

q=l 9=1 

Assuming that Cq has a mild g-dependence, this expansion converges when the "effective" 
expansion parameter gM/r+^^T in(2K) ^ i jg small. Therefore, the convergence region is given 
by 

/i/T < -iVTln(2K). (14) 

Since actual simulations are performed at k > 0.125 (ln(2K) < —1.386 ■ ■ ■ )' ^^e convergence 
region is much wider than that of Taylor expansion, fi/T < 1. This analysis also indicates 
that the convergence becomes better for larger iVr for an fixed. 



4 



C. An upper bound for the absolute value of the phase 



The origin of the sign problem is a large fluctuation of the phase. Understanding the 
properties of the phase in more detail may open a new insight on either solving or avoiding 
the sign problem. Here we derive an analytical upper bound for the leading term of the phase 
in the hopping parameter expansion. While this is just an approximated upper bound, it 
reveals interesting features of the phase, which we later confirm by numerical simulations. 

To derive the upper bound, we use the inequality + < + for x,y G C repeatedly 
to find that 



1^1 < ^2sinh(g/x/T)|Im[\/(9)]| < ^ 2 sinh(g/i/T)|V('')|. 



(15) 



q=l 



To obtain an explicit lattice parameter dependence, let us apply some approximations. First 
we truncate the winding number expansion at g = 1. In the previous section, this truncation 
has been shown to be valid for small k, and large N^. After this truncation, the upper bound 
becomes 



J] 2 sinh (g/i/T) I V^' 



> 2sinh(/i/r)|y(^)|. (16) 
We then use the explicit expressions for 1/(1) given in eqs. flA33|) and (IA34|) . and obtain 



|\/«| = |tr(i7+)|+0((2/€)2^T) 
= (2«:)^--trspin[(p(^))^-] 
< (2k)^t ■ 2 ■ 3Nl, 



^ ^ trpolor 



_X4=1 



(17) 



where in the last line we used the fact that the trace of an SU(3) matrix satisfies |trcoiorf^| 



approximate upper bound given by 



3 and trspin-P*-^^ 



2. Assembling various pieces, we obtain an 



1^1 < 12Nl{2K)^^smh{fi/T). 



(18) 



Consider the case that we keep fi/T fixed. Then this expression tells us that the phase 
becomes exponentially suppressed for larger temporal extent Nt, while it increases in pro- 
portion to the lattice spatial volume A'^l- We also observe that the phase vanishes in the 
static limit k — )■ with spatial volume fixed. Therefore the order of the static limit and the 
thermodynamic limit is subtle. 

Alternatively, one may wish to consider the case that a/i is kept fixed. In this case, we 
rewrite /i/T = afiN^ in eq. f ll8p and find that 



- 6NI{2k)^'^ expiafiN^) 
= 6Nlexp{{ln{2K) + afijN 



T 



(19) 



The bound still decreases exponentially for large A^x if < — ln(2K) is satisfied. This 
bound on a/x is mild since — ln(2fi:) > 1.386 ■ ■ ■ for n > 0.125, and hence covers the region 
/i~l/a~l — 2 GeV in current typical simulations. Suppose that we change the temperature 
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as a function of A^t while keeping other parameters a, A'l, k, and a/x fixed. By increasing A^t 
the upper bound in eq. (fTO!) becomes exponentially suppressed and the sign problem may 
be avoided at low temperatures. 

We have seen that the magnitude of the phase becomes smaller as the temporal lattice 
size A^T increases. We still have to examine whether the magnitude of the phase is under 
control even if the spatial volume becomes large. To this end, let us consider the situation 
that all physical scales, viz., 

• T = l/(aA'T): temperature, 

• yU.: chemical potential, 

• V = Li^ = [aNi^Y: spatial volume, 

• m: physical mass. 

are kept fixed. If we change the lattice spacing by a factor 6 > 0, then we have to rescale 
the other parameters to keep all the physical scales constant: 

• a — )■ a/h, 

• Nt — )■ BNt, 

• iVL hNi^, 

• a/i — )■ afi/b, 

• K, — y K,' ^ K, [am — > am/b), 

where we have made use of the fact that k is not so sensitive to the change of the lattice 
spacing compared to other parameters. The ratio of the bound before and after the scaling 
is given by 



bound of 


1^ 


1 after 


bound of 


^1 


before 



12{bNi,f{2K')''^^ sinh(/i/T) 
12iV3(2K)^Tsinh(/i/T) 



6^(2/€)^T{b-i)_ (20) 



For typical values of 2k and A'^x (say 2k, < 1 and A^t > 4), reducing the lattice spacing also 
reduces the magnitude of the phase since the exponential factor (2k)^t(^~^) dominates over 
the polynomial factor b^. It is important to check whether this behavior of the phase is 
observed for the exact phase. This is under investigation with numerical simulations. 



III. NUMERICAL TEST IN iV/ = 4 QCD 
A. Simulation parameters 

We employ the clover- improved Wilson quark action for four degenerate flavors and the 
Iwasaki gauge action. The simulation parameters are the same as in Ref. [li, namely 



/3 = 1.60, K = 0.1371 and = 1.9655, which corresponds to the lattice spacing a = 0.328 
fm and m^r = 834 MeV. The spatial lattice size is set to A'l = 6 which gives L = Ni^a ~ 2 
fm, while we employ various temporal sizes of A'r = 4, 6, 8 and 12 which correspond to the 
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temperature T = l/(aA'T) = 50 — 150 MeV. The quark chemical potential is also varied in 
the range 0.05 < a/i < 0.8. 

We use the conventional HMC algorithm for the phase- quenched quark determinant with 
the iso-spin chemical potential /iu = —/id- Two independent pseudo-fermions are prepared 
to incorporate A'^f = 4 dynamical quarks. We set the trajectory length to unity and vary the 
step size dr = 1/50 — 1/240 depending on a/i and A^t such that the HMC acceptance rate 
stays around 90 %. Since the fermion force tends to be large for larger chemical potential, 
we diminish dr to keep the HMC acceptance rate. For each parameter set, 5000 — 20000 
trajectories are accumulated. We employ the jackknife analysis for the error estimate. 



B. Convergence check and analysis of the exact phase 

We first check the exponential decay behavior of V^'^^ as a function of the winding index q 
expected from eq. f|T2|) . Figure [1] shows the q and dependence of the real and imaginary 
parts of V^'^'^ at a/i = 0.2 measured on 200 phase-quenched configurations. We observe a 
clear exponential decay behavior, though the decay rate is somewhat milder than eq. ( !T2|) . 

We next examine the behavior of the phase. In Fig. [2] we plot the q dependence of 
2 sinh(g/i/T)|Im[y^'')]| in eq. fll5p . Even after the factor sinh(g/i/T) is combined, the results 
show exponential fall-off as a function of q for iVr > 6. As expected, the convergence for 
the phase becomes better for larger A^x- 

Finally let us examine the exact phase by calculating it directly from the quark deter- 
minant. The calculational procedure is explained in Subsection IIII CI together with that of 
quark number density. Figure [3] show the distribution of the phase for the single Wilson- 
Dirac determinant measured on the same phase-quenched configurations as in Figs. [1] and 
[21 We find that the shape of the distribution becomes sharper for larger A'r. 

We fit the distribution of phase by a normal distribution and extract the width a as a 
function of a/i and N^^, which is plotted in Fig. HI The exponential decrease of the width 
nicely confirms our expectation. The slope of the exponential decrease is similar to that 
calculated from the upper bound estimate in eq. (IT^ . 

Another important quantity is the averaged phase factor (e'^^)|| = (cos(4^))|| where the 
factor 4 accounts for the degenerate A^f = 4 flavors. The average is taken on the phase- 
quenched configurations and the precise definition of (...)|| is again given in Subsection IIII C[ 
Figure [5] shows the result as a function of a/i for A^t = 4, 6, 8, 12 for A^'l = 6. The average 
stays closer to unity for larger A^t, even for rather large a/i. Note that the charged pion 



condensation [13| could occur when a/i > a/ic = ^ 0.7 at T = 0. 



C. Physical observables with reweighting method 

Since the phase can be controlled by increasing the temporal size of the lattice, we 
expect that the phase reweighting method would work in the low temperature region. The 
reweighting formula for the degenerate A^f-fiavor case is given by 

iO) - (21) 
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where the reweighted and the phase- quenched ensembles averages are defined as 



^ J[dU]e-'o{detDr^O[U] 
^ ' J[dU]e-3o(detD)^^ ' ^ ^ 

J[dU]e-'o\detD\^^0[U] 



We measure the quark number, the plaquette, the polyakov loop, and their susceptibilities 
by using the reweighting formula in eq. ( 12T|) . The plaquette and the polyakov loop are 
measured at every trajectory, while the quark number nq and the susceptibility Xq ^ire 
exactly calculated for every 10 trajectories without relying on the noise method. 

The expressions for riq and Xq are derived from the quark determinant in the reduced 
form given by 

det D{fi) = Ao det[l - Ho - e^/^H^ - e^^/^H^], (24) 
where with k = 0,± are given in eqs. flA25p — flA27p . Differentiating the grand canonical 



partition function with respect to /i/T, we derive 

^ = ^(t.iA'-'A1). (25) 

N'^ / / . \ 2\ /V^ / . \ 2 



with 



K{fi/T) = l-Ho-e^'/^H+-e~''/'^H_, (27) 
kifi/T) = -[e'^/^if+-e-''/^iJ„], (28) 
K{fi/T) = -[e^/^if+ + e-^/^if_]. (29) 



The phase is given by 



Re [det K] ^ ' 

Let us make a comment on the computational cost. The matrices H±, Hq are dense 
matrices of rank 12Nf^. Once they are calculated and stored in memory, the phase and the 
quark number and its susceptibility can be calculated exactly up to machine precision. The 
computational cost of Hk, in turn, scales as (7Vl(-^t/2 — 1))'^ due to the LU decomposition 
for the inversion of the matrices and D(3) in eq. (jATOj) which have rank 12A'l(A'^t/2 — 1). 
It is crucial to reduce these inversion costs, and we have implemented a further reduction, 
achieving a computational cost proportional to NlN^. The required memory scales as N^. 

Figure El shows our results for the quark number and its susceptibility for A'^x = 4. 
Reweighted averages are shown by red symbols and connected by red solid line ("with 
phase"), whereas the phase quenched averages are shown in green symbols and connected 
by green dashed line ("without phase"). We observe that the effect of the phase is not very 
large for a/i up to 0.6 over which the average of the phase stays non-zero (see Fig. [5]). The 
peak of the susceptibilities slightly shift to lower a/x for the reweighted average. A similar 
phenomenon is observed also for the plaquette and the polyakov loop shown in Figs. [7] and 
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[HI respectively. At a/x = 0.7 and beyond, the normal average loses control since the phase 
average becomes consistent with zero. 

The presence of the peak observed in the susceptibilities suggests a phase transition. 
Numerically the peaks are located at a/i ~ 0.21 where n^/T^ ~ 7.7. In terms of the 
variables Hb/T, where /ib = 3/i is the baryon chemical potential, and the baryon number 



ub used in Fig. 7 (T = 0.95Tc) of Ref. [12|, these values translate into Hb/T ~ 2.5 and 



ub ~ 8.6, which are roughly consistent with the region of S shape found by the canonical 
approach in Ref. [ll ]. 

The histograms of the plaquette on the phase-quenched configuration are shown in Fig. |9l 
There is a double peak around a/i = 0.21 indicating a first order phase transition. In order 
to pin down the order of transition, one has to carry out a finite size scaling. This is left for 
future work. 

We have also carried out the reweighting simulations increasing the temporal size from 
iVx = 4 to 6, 8 and 12, thus lowering the temperature from T = 150 MeV to lOOMeV down 
to 50 MeV. We could calculate observables reliably over a significant range of a/i since the 
phase average stays non-vanishing as shown in Fig. |3 The results, however, showed only 
smooth variations, and we did not observe signals of the phase transition over the region 
of afi where the phase is under control. It may well be that, at lower temperatures, the 
expected phase transition takes place at larger afi where one runs into the sign problem, 
i.e., the phase average becomes very small or even vanishes. This possibility reminds us of 
Ref. jsj in which the authors showed for the two-fiavor random matrix model that the phase 
average vanishes at the boundary of a region surrounding the phase transition line in the 
(T, /i) plane. 



IV. CONCLUDING REMARKS 

We studied the phase of the quark determinant in QCD for finite chemical potential. 
Analytical estimates show that the winding number expansion has a nice convergence prop- 
erty, and we have confirmed this through a numerical calculation of the expansion terms 
in phase-quenched QCD. We hope that this expansion proves to be useful to explore the 
high density region at low temperatures, which has not been explored well by numerical 
simulations. 

We also found that the magnitude of the phase becomes smaller for larger temporal 
size of the lattice. We have shown this both by analytical estimates based on the hopping 
parameter expansion and by an exact evaluation of the determinat. The quark mass which we 
employ in this paper is quite heavy, however. Hence we have to check whether this property 
persists toward the physical quark mass region. As a further result, simple scaling argument 
combined with the hopping parameter expansion of the phase leads to the prediction that 
the phase becomes smaller for smaller lattice spacing with all the physical scales fixed. This 
is an interesting possibility which we hope to check in future. 

Finally, we explored the phase structure for A'^f = 4 QCD with the use of the reweighting 
method, and found signals indicative of a first order phase transition on an 6^ x 4 lattice 



at T ^ 150 MeV and /i !^ 110 MeV. This result is consistent with that of Ref. [12| using 
canonical methods. Needless to say, in order to pin down the order of the transition, one 
has to carry out the finite size scaling. 
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Appendix A: Reduction for the quark determinant of the Wilson-Dirac operator 



To fix and summarize our notations, we briefly review a reduction technique in Ref. [11 



1. Wilson-Dirac operator 

We consider the Wilson-Dirac operator^ on the lattice with the size A^l ^ -^t- The Wilson 
fermion action is given by 

S'w = ^ (^(x)^(x) 

X 

3 

-2k ^[^(x)P!^^t/(x, k)^{x + k)+ ^{x + k)Pi^^U{x, A:) 

k=l 

-2K[e"^^/^(x)P!^^f/(x, 4)iP{x + 4) + e-''^'ip{x + 4)P|'^^f/(x, 4)1"V^(x)]) , (Al) 



with the projections P^^ = (l±7i,)/2. The fermion fields satisfy the anti-periodic boundary 
conditions 

V^(x,A^T + l) = -^(x,l), (A2) 

V^(x,A^T + l) = -^(x,l). (A3) 

Applying the follwing transformation to the fermion fields, 

V^(x,a;4) — > e-"^^^7^(x,X4), (A4) 

z^(x,X4) ^ e'*^^4^A(x,X4), (A5) 

the yU-dependence is localized in the time-hopping term at X4 = 1 and X4 = At planes. 
By decomposing the above lattice in the time direction as illustrated in Fig. [TOl 

1,2,3,..,A^t/2-1, (A6) 

Nt/2, (A7) 

A^t/2 + 1,A^t/2 + 2,...,A^t- 1, (A8) 

A^T, (A9) 



A(i) 


for 


X4 


A(2) 


for 


X4 


A(3) 


for 


X4 


A(4) 


for 





^ Inclusion of the clover term is straightforward. 
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the decomposed operator can be written as (time blocked form) 



^(1) 


^(12) 







^(21) 


Dm 


^(23) 








D{32) 




-0(34) 







^(43) 


^(4) 



(AlO) 



where the //-dependence is shown exphcitly. Note that and -D(jj) do not depend on n 
and follow the 75-hermiticity and 75-relation 



for i = 1,2,3,4, 



75^(i)75, 

75D(,,)75, for (tj) = (12), (23), (34), (41). 



(All) 
(A12) 



This decomposition is useful for the factorization of fermion determinant as we will see later. 



2. Reduction for the time direction 

By making use of the following formula repeatedly 



det 



A B 
C D 



det [A] det[D] det[l - D-^CA-^B], 



(A13) 



one can obtain the reduced expression of the determinant 
det -D(/i) = det D(i) det D(3) det -D(2*2) det -D(4*4) 

(4*4) 

where we have defined 



X detil - D^.L^(e'^/^D(412) + D(432))/^("2l2)(e"^/^^(214) 



D 
D 
D 
D 
D 
D 



(2*2) 
(4*4) 
(412) 
(432) 
(214) 
(234) 



^(2) - /^(21)^(i5^(12) - 15(23)^(3)^(32) 



D 



(4) 



(1) 

^(41)^(^l!^(14) 



(3)^ 

D ^43) D^s]D (34), 



^(41)^(^1)^(12), 
^(43) ^("s! -0(32), 

^(21)^(^15-0(14), 
-0(23) -0(^35 -0(34)- 



D(234))], (A14) 
(A15) 

(A16) 
(A17) 
(A18) 
(A19) 
(A20) 



From eqs. (lAlip and (1A12P they have the 75-hermiticity and the 75-relation, 

75-D(i*i)75 for z = 2,4, 



(i*t) 



'l5D(jki)'j5 for i,j = 2, 4 and k = l,3. 
The expression of eq. (lA14p is simplified as 

det D(/i) = Aodet[l- Ho- e^'^H^ - e-^^'^H^ 



(A21) 
(A22) 

(A23) 
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with the introduction of 



Aq = det _D(i) det D(s) det -D(2*2) det -D(4*4), (A24) 

^0 = ^(4l4)^{412)^(2l2)^{214) + ^(4l4)^{432)^(2l2)^(234), (A25) 

H+ = Z^(4i,)Z^(412)/^(2l2)^{234), (A26) 

H~ = -0(4^4) -D(432) £'(2^2) -^(214)- (A27) 



From eqs. f lAll|) and flA21|) . Aq is shown to be real. An geometrical meaning of Hk {k = 0, ±) 



is as follows (see also Fig. [TT]) : 

• Hq: paths through domains A(4) — )■ A(i)(orA(3)) — A(2) — ?■ A(i)(orA(3)) — A(4) 
=^ NO winding around time direction. 

• H^: paths through domains A(4) — )■ A(i) A(2) — > A(3) A(4) 
=^ Forward winding around time direction. 

• H^: paths through domains A(4) — )■ A(3) A(2) — t- A(i) — j- A(4) 

Backward winding around time direction. 



From eqs. (]A21|) and ( ]A22p we find that they have 75-relation 

Hi = 75l^(4*4)i^-.I^(-4l4)75 (A28) 

for k = 0,±. 

3. Winding number expansion 

We can rewrite det -D(/i) in eq. flA23p as follows, 

detDifi) = Aq det[l - Hq - e^"'^ H+ - e-^^'^H^] (A29) 
= Aq exp (Tr ln[l - Hq - e''/^H+ - e-^^'^HJ]) (A30) 

= Ao exp ( - V -Tr[ifo + e^/^if+ + e-^'/^Hr \ (A31) 

V / 

(00 00 ^ \ 

^ Tt[Hu,Hu,...H,U , (A32) 
ggZ n=l fcl+...+fc„=(J / 

where ki G {0,±1} and g is an integer which counts how many times the individual loop 
Tr[iffcjif/c2...iffc^] winds around the time direction. In order to simplify the equation, we 
introduce 

00 

yi,) = Y,V^'\ (A33) 

n=l 

V^'^ = - E ^4Hk,H,,...Ha (A34) 



n 

ki+...+kn=q 
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where ki G {0,±1}. Note that Vn'^^ = for q > n by definition. Since for an actual 
evaluation of V^'^'^ one cannot of course sum up all the terms in eq. flA33l) . we adopt the 



truncation scheme in Ref. 11 



V^(°) = tT[Ho] + ^tT[{Hof]+tT[H+H. 



for g = 1, 2, 3, ... 



(A35) 
(A36) 



The relative truncation errors are 0((2k)^'^) thus can be safely neglected. From eg. ( 1A28I1 
one find 

Ti[Hk,Hk,...HkJ = Tr[i/„fc„...iJ„fc,ff„fcJ*. (A37) 

This yields 



yi-i) — yii)* _ 



(A38) 
(A39) 



Thus one can show eqs. ([2]) and ([3]). Finally we summarize the winding number expansion 
of the logarithm of the determinant 



det D{fi) = Aoexp 



(A40) 
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Nt=4 



Nt=6 



Nt=8 



Nt=12 




FIG. 1: Absolute values of real (blue points) and imaginary (red points) parts of V^'^^ for 
q = 1,2,..., 10 with various A^t- They are measured at afi = 0.2 employing 200 phase-quenched 
configurations. Red line denotes an expected asymptotic behavior (2k)'?^t given in eq. (fT^ . 



Nt=4 Nt=6 Nj=8 Nt=12 




FIG. 2: Absolute value of each term of the phase in eq. dS]): 2sinh(g/x/T)|Im[l/*^''-']| foiq = 1,2, ...,10 
with various A^t on the same configurations as those in Fig. [TJ 
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FIG. 3: Histogram of the phase at afi = 0.2 with various Nj^. Other parameters are the same as 
those of Fig. [TJ 




FIG. 4: Width of the distribution of phase obtained by fitting the distribution to the normal 
distribution as a function of the temporal size A't for afi = 0.1 to 0.5. Other parameters are the 
same as those of Fig. [TJ 
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FIG. 5: Averaged reweighting factor as as function of afj, with various Nt^. Other parameters are 
the same as those in Fig. [TJ 
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FIG. 6: Quark number density (upper) and its susceptibility (lower) as a function of a// (or the 
baryon chemical potential = 3/u in unit of the temperature for upper horizontal axis) with and 
without the phase factor. Light blue dotted arrow shows the location of transition point determined 
by the canonical approach [1^. Red solid arrow shows the location of a peak of the quark number 
susceptibility (a/i = 0.21). The resluts at afi = 0.7 and 0.8 have the huge error bar because of the 
extremely small reweighting factor (an appearance of the sign problem) as shown in Fig. [SJ 
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FIG. 7: Averaged plaquette value (upper) and its susceptibility (lower) as a function of a/x with 
and without phase factor, namely quark and iso-spin chemical potential cases. 
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FIG. 8: Averaged real part of the polyakov loop (upper) and its susceptibility (lower) as a function 
of a/i with and without the phase factor. 
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FIG. 9: Histogram of plaquette value on the phase-quenched configuration for A^t = 4 with various 
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FIG. 10: Domain decomposition in two dimensional case. Time index X4 runs from 1 to A^t- 
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FIG. 11: Geometrical pictures of Hq, H±. 
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